{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018  by D. Koehn, heterogeneous models are from [this Jupyter notebook](https://nbviewer.jupyter.org/github/krischer/seismo_live/blob/master/notebooks/Computational%20Seismology/The Finite-Difference Method/fd_ac2d_heterogeneous.ipynb) by Heiner Igel ([@heinerigel](https://github.com/heinerigel)), Florian Wölfl and Lion Krischer ([@krischer](https://github.com/krischer)) which is a supplemenatry material to the book [Computational Seismology: A Practical Introduction](http://www.computational-seismology.org/), notebook style sheet by L.A. Barba, N.C. Clementi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href=\"https://fonts.googleapis.com/css?family=Merriweather:300,300i,400,400i,700,700i,900,900i\" rel='stylesheet' >\n",
       "<link href=\"https://fonts.googleapis.com/css?family=Source+Sans+Pro:300,300i,400,400i,700,700i\" rel='stylesheet' >\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' >\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: rgb(245,245,245);\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width */\n",
       "    width: 800px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    width: 1000px;\n",
       "    margin: auto;\n",
       "    padding-left: 0em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "margin-top:0.5em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    border: 1px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\n",
       "div.cell.code_cell {\n",
       "    background-color: rgb(256,256,256); \n",
       "    border-radius: 0px; \n",
       "    padding: 0.5em;\n",
       "    margin-left:1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    line-height: 140%;\n",
       "    font-size: 110%;\n",
       "    width:680px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-style:regular;\n",
       "    font-weight: bold;    \n",
       "    font-size: 250%;\n",
       "    line-height: 100%;\n",
       "    color: #004065;\n",
       "    margin-bottom: 1em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: bold; \n",
       "    font-size: 180%;\n",
       "    line-height: 100%;\n",
       "    color: #0096d6;\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "\tfont-size: 150%;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "    color: #008367;\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: 300; \n",
       "    font-size: 100%;\n",
       "    line-height: 120%;\n",
       "    text-align: left;\n",
       "    width:500px;\n",
       "    margin-top: 1em;\n",
       "    margin-bottom: 2em;\n",
       "    margin-left: 80pt;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    font-weight: regular;\n",
       "    font-size: 130%;\n",
       "    color: #e31937;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .5em;\n",
       "    margin-top: 1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'Source Code Pro', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 9pt;\n",
       "    line-height: 100%;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "\t\n",
       "    \n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"], \n",
       "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise: 2D acoustic FD modelling of the Marmousi-2 model\n",
    "\n",
    "In this exercise, you have to apply all the knowledge about FD modelling, we covered so far. While the modelling examples in the last lesson where quite simple, we now calculate the 2D acoustic wave propagation in a more realistic problem called the **Marmousi-2 model**.\n",
    "Developed in the 1990s by the French Petroleum Institute (IFP) ([Versteeg, 1994](https://library.seg.org/doi/abs/10.1190/1.1437051)), the Marmousi model is a widely used benchmark problem for seismic imaging and inversion techniques. Beside the original acoustic version of the model, an elastic version was developed by [Martin et al. (2006)](https://library.seg.org/doi/abs/10.1190/1.2172306).\n",
    "\n",
    "The Marmousi-2 model consists of a 460 m thick water layer above an elastic subseafloor model.  The sediment model is very simple near the left and right boundaries but rather complex in the centre. At both sides, the subseafloor is approximately horizontally layered, while steep thrust faults are disturbing the layers in the centre of the model. Embedded in the thrust fault system and layers are small scale hydrocarbon reservoirs.\n",
    "\n",
    "##### Exercise\n",
    "\n",
    "Setup and model the 2D acoustic wave propagation in the Marmousi-2 model:\n",
    "\n",
    "- Define the model discretization based on the given Marmousi-2 P-wave velocity model. Maximum wave propagation time should be 6 s \n",
    "- Calculate the central frequency $f_0$ of the source wavelet based on the grid dispersion criterion \n",
    "\n",
    "\\begin{equation}\n",
    "dx \\le \\frac{vp_{min}}{N_\\lambda f_0}, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "which you can use for the FD modelling run, based on the pre-defined $dx$ of the Marmousi-2 model , minimum P-wave velocity $vp_{min}$ and $N_\\lambda = 4$ gridpoints per dominant wavelength.\n",
    "- Start a modelling run for the Marmousi-2 model with the 3-point spatial FD operator. Place an airgun for a central shot at x = 5000 m at a depth = 40 m below the sea surface. Do not forget to calculate an appropriate time step $dt$\n",
    "- Add an additional function `update_d2px_d2pz_5pt` to approximate the 2nd spatial derivatives by the 5-point FD operator in the 2D acoustic FD code derived [here](http://nbviewer.jupyter.org/github/daniel-koehn/Theory-of-seismic-waves-II/blob/master/04_FD_stability_dispersion/3_fd_taylor_operators.ipynb). Add an option `op` to switch between the 3-point and 5-point operator\n",
    "- Start an additional modelling run for the Marmousi-2 model with the 5-point operator.\n",
    "- Imagine you place an Ocean-Bottom-Cable (OBC) on the seafloor of the Marmousi-2 model. Calculate an OBC shot gather, by placing receivers at each gridpoint of the Cartesian model in x-direction in a depth of 460 m. Modify the FD code to record seismograms at each receiver position. Do not forget to return the seismograms from the modelling function.\n",
    "- Plot and compare the seismograms produced by the 3- and 5-point operators, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Import Libraries \n",
    "# ----------------\n",
    "import numpy as np\n",
    "from numba import jit\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from pylab import rcParams\n",
    "\n",
    "# Ignore Warning Messages\n",
    "# -----------------------\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The elastic version of the Marmousi-2 model, together with FD modelled high frequency shot gather data is available from [here](http://www.agl.uh.edu/downloads/downloads.htm). The central part of the P-wave velocity model of the Marmousi-2 with the spatial discretization: \n",
    "\n",
    "$nx = 500$ gridpoints\n",
    "\n",
    "$nz = 174$ gridpoints\n",
    "\n",
    "$dx = dz = 20\\; m$\n",
    "\n",
    "is available as IEEE little endian binary file `marmousi_II_marine.vp` in the `marmousi-2` directory. It can be imported to Python with the following code snippet:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Marmousi-2 Vp model\n",
    "# --------------------------\n",
    "\n",
    "# DEFINE MODEL DISCRETIZATION HERE!\n",
    "nx =       # number of grid points in x-direction\n",
    "nz =       # number of grid points in z-direction\n",
    "dx =       # spatial grid point distance in x-direction (m)\n",
    "dz = dx        # spatial grid point distance in z-direction (m)\n",
    "\n",
    "# Define model filename\n",
    "name_vp = \"marmousi-2/marmousi_II_marine.vp\"\n",
    "\n",
    "# Open file and write binary data to vp\n",
    "f = open(name_vp)\n",
    "data_type = np.dtype ('float32').newbyteorder ('<')\n",
    "vp = np.fromfile (f, dtype=data_type)\n",
    "\n",
    "# Reshape (1 x nx*nz) vector to (nx x nz) matrix \n",
    "vp = vp.reshape(nx,nz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After reading the model into Python, we can take a look at it ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot Marmousi-2 vp-model\n",
    "# ------------------------\n",
    "\n",
    "# Define xmax, zmax and model extension\n",
    "xmax = nx * dx\n",
    "zmax = nz * dz\n",
    "extent = [0, xmax, zmax, 0]\n",
    "\n",
    "fig = plt.figure(figsize=(12,3))  # define figure size\n",
    "\n",
    "image = plt.imshow((vp.T)/1000, cmap=plt.cm.viridis, interpolation='nearest', \n",
    "                   extent=extent)\n",
    "\n",
    "cbar = plt.colorbar(aspect=10, pad=0.02)\n",
    "cbar.set_label('Vp [km/s]', labelpad=10)\n",
    "plt.title('Marmousi-2 model')\n",
    "plt.xlabel('x [m]')\n",
    "plt.ylabel('z [m]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [],
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Definition of modelling parameters\n",
    "# ----------------------------------\n",
    "# DEFINE MAXIMUM RECORDING TIME HERE!\n",
    "tmax =  # maximum wave propagation time (s)\n",
    "\n",
    "# DEFINE YOUR SHOT POSITION HERE!\n",
    "xsrc =   # x-source position (m)\n",
    "zsrc =   # z-source position (m)\n",
    "\n",
    "# CALCULATE DOMINANT FREQUENCY OF THE SOURCE WAVELET HERE!\n",
    "f0   = # dominant frequency of the source (Hz)\n",
    "print(\"f0 = \", f0, \" Hz\")\n",
    "t0   = 4.0/f0   # source time shift (s)\n",
    "\n",
    "isnap = 2  # snapshot interval (timesteps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit(nopython=True) # use JIT for C-performance\n",
    "def update_d2px_d2pz_3pt(p, dx, dz, nx, nz, d2px, d2pz):\n",
    "    \n",
    "    for i in range(1, nx - 1):\n",
    "        for j in range(1, nz - 1):\n",
    "                \n",
    "            d2px[i,j] = (p[i + 1,j] - 2 * p[i,j] + p[i - 1,j]) / dx**2                \n",
    "            d2pz[i,j] = (p[i,j + 1] - 2 * p[i,j] + p[i,j - 1]) / dz**2\n",
    "        \n",
    "    return d2px, d2pz "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ADD THE SPATIAL 5-POINT OPERATORS HERE!\n",
    "@jit(nopython=True) # use JIT for C-performance\n",
    "def update_d2px_d2pz_5pt(p, dx, dz, nx, nz, d2px, d2pz):\n",
    "    \n",
    "        \n",
    "    return d2px, d2pz "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define simple absorbing boundary frame based on wavefield damping \n",
    "# according to Cerjan et al., 1985, Geophysics, 50, 705-708\n",
    "def absorb(nx,nz):\n",
    "\n",
    "    FW = 60     # thickness of absorbing frame (gridpoints)    \n",
    "    a = 0.0053\n",
    "    \n",
    "    coeff = np.zeros(FW)\n",
    "    \n",
    "    # define coefficients in absorbing frame\n",
    "    for i in range(FW):    \n",
    "        coeff[i] = np.exp(-(a**2 * (FW-i)**2))\n",
    "\n",
    "    # initialize array of absorbing coefficients\n",
    "    absorb_coeff = np.ones((nx,nz))\n",
    "\n",
    "    # compute coefficients for left grid boundaries (x-direction)\n",
    "    zb=0 \n",
    "    for i in range(FW):\n",
    "        ze = nz - i - 1\n",
    "        for j in range(zb,ze):\n",
    "            absorb_coeff[i,j] = coeff[i]\n",
    "\n",
    "    # compute coefficients for right grid boundaries (x-direction)        \n",
    "    zb=0\n",
    "    for i in range(FW):\n",
    "        ii = nx - i - 1\n",
    "        ze = nz - i - 1\n",
    "        for j in range(zb,ze):\n",
    "            absorb_coeff[ii,j] = coeff[i]\n",
    "\n",
    "    # compute coefficients for bottom grid boundaries (z-direction)        \n",
    "    xb=0 \n",
    "    for j in range(FW):\n",
    "        jj = nz - j - 1\n",
    "        xb = j\n",
    "        xe = nx - j\n",
    "        for i in range(xb,xe):\n",
    "            absorb_coeff[i,jj] = coeff[j]\n",
    "\n",
    "    return absorb_coeff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FD_2D_acoustic code with JIT optimization\n",
    "# -----------------------------------------\n",
    "def FD_2D_acoustic_JIT(vp, dt,dx,dz,f0,xsrc,zsrc,op):        \n",
    "    \n",
    "    # calculate number of time steps nt \n",
    "    # ---------------------------------\n",
    "    nt = (int)(tmax/dt)\n",
    "    \n",
    "    # locate source on Cartesian FD grid\n",
    "    # ----------------------------------\n",
    "    isrc = (int)(xsrc/dx)  # source location in grid in x-direction\n",
    "    jsrc = (int)(zsrc/dz)  # source location in grid in x-direction    \n",
    "    \n",
    "    # Source time function (Gaussian)\n",
    "    # -------------------------------\n",
    "    src  = np.zeros(nt + 1)\n",
    "    time = np.linspace(0 * dt, nt * dt, nt)\n",
    "\n",
    "    # 1st derivative of Gaussian\n",
    "    # src  = -2. * (time - t0) * (f0 ** 2) * (np.exp(- (f0 ** 2) * (time - t0) ** 2))\n",
    "    \n",
    "    # 2nd derivative of a Gaussian\n",
    "    src = (1.0 - 2.0*(np.pi**2)*(f0**2)*((time-t0)**2)) * np.exp(-(np.pi**2)*(f0**2)*((time-t0)**2))\n",
    "    \n",
    "    # define clip value: 0.1 * absolute maximum value of source wavelet\n",
    "    clip = 0.5 * max([np.abs(src.min()), np.abs(src.max())]) / (dx*dz) * dt**2\n",
    "    \n",
    "    # Define absorbing boundary frame\n",
    "    # -------------------------------    \n",
    "    absorb_coeff = absorb(nx,nz)\n",
    "    \n",
    "    # Define squared vp-model\n",
    "    # -----------------------        \n",
    "    vp2 = vp**2\n",
    "    \n",
    "    # Initialize empty pressure arrays\n",
    "    # --------------------------------\n",
    "    p    = np.zeros((nx,nz)) # p at time n (now)\n",
    "    pold = np.zeros((nx,nz)) # p at time n-1 (past)\n",
    "    pnew = np.zeros((nx,nz)) # p at time n+1 (present)\n",
    "    d2px = np.zeros((nx,nz)) # 2nd spatial x-derivative of p\n",
    "    d2pz = np.zeros((nx,nz)) # 2nd spatial z-derivative of p \n",
    "    \n",
    "    # INITIALIZE SEISMOGRAMS HERE! \n",
    "    # ----------------------------\n",
    "        \n",
    "    # Initalize animation of pressure wavefield \n",
    "    # -----------------------------------------    \n",
    "    fig = plt.figure(figsize=(7,3))  # define figure size\n",
    "    extent = [0.0,xmax,zmax,0.0]     # define model extension\n",
    "    \n",
    "    # Plot Vp-model\n",
    "    image = plt.imshow((vp.T)/1000, cmap=plt.cm.gray, interpolation='nearest', \n",
    "                        extent=extent)    \n",
    "    \n",
    "    # Plot pressure wavefield movie\n",
    "    image1 = plt.imshow(p.T, animated=True, cmap=\"RdBu\", alpha=.75, extent=extent, \n",
    "                          interpolation='nearest', vmin=-clip, vmax=clip)    \n",
    "    plt.title('Pressure wavefield')\n",
    "    plt.xlabel('x [m]')\n",
    "    plt.ylabel('z [m]')\n",
    "           \n",
    "    plt.ion()    \n",
    "    plt.show(block=False)\n",
    "    \n",
    "    # Calculate Partial Derivatives\n",
    "    # -----------------------------\n",
    "    for it in range(nt):\n",
    "    \n",
    "        # FD approximation of spatial derivative by 3 point operator\n",
    "        if(op==3):\n",
    "            d2px, d2pz = update_d2px_d2pz_3pt(p, dx, dz, nx, nz, d2px, d2pz)\n",
    "        \n",
    "        # ADD FD APPROXIMATION OF SPATIAL DERIVATIVES BY 5 POINT OPERATOR HERE!\n",
    "        #if(op==5):\n",
    "\n",
    "        # Time Extrapolation\n",
    "        # ------------------\n",
    "        pnew = 2 * p - pold + vp2 * dt**2 * (d2px + d2pz)\n",
    "\n",
    "        # Add Source Term at isrc\n",
    "        # -----------------------\n",
    "        # Absolute pressure w.r.t analytical solution\n",
    "        pnew[isrc,jsrc] = pnew[isrc,jsrc] + src[it] / (dx * dz) * dt ** 2\n",
    "        \n",
    "        # Apply absorbing boundary frame\n",
    "        # ------------------------------\n",
    "        p *= absorb_coeff\n",
    "        pnew *= absorb_coeff\n",
    "        \n",
    "        # Remap Time Levels\n",
    "        # -----------------\n",
    "        pold, p = p, pnew\n",
    "        \n",
    "        # WRITE SEISMOGRAMS HERE!\n",
    "    \n",
    "        # display pressure snapshots \n",
    "        if (it % isnap) == 0:            \n",
    "            image1.set_data(p.T)\n",
    "            fig.canvas.draw()\n",
    "    \n",
    "    # DO NOT FORGET TO RETURN THE SEISMOGRAM HERE!      "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model wave propagation in the Marmousi-2 model with 2D acoustic FD code\n",
    "\n",
    "Time to model acoustic wave propagation in the Marmouisi-2 model using the 3-point operator. We only have to define the timestep $dt$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run 2D acoustic FD modelling with 3-point spatial operater\n",
    "# ----------------------------------------------------------\n",
    "%matplotlib notebook\n",
    "op = 3  # define spatial FD operator (3-point) \n",
    "\n",
    "# DEFINE TIME STEP HERE!\n",
    "dt = # time step (s)\n",
    "\n",
    "FD_2D_acoustic_JIT(vp,dt,dx,dz,f0,xsrc,zsrc,op)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run 2D acoustic FD modelling with 5-point spatial operater\n",
    "# ----------------------------------------------------------\n",
    "%matplotlib notebook\n",
    "op = 5  # define spatial FD operator (5-point) \n",
    "\n",
    "# DEFINE TIME STEP HERE!\n",
    "dt   =  # time step (s)\n",
    "\n",
    "FD_2D_acoustic_JIT(vp,dt,dx,dz,f0,xsrc,zsrc,op)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "# PLOT YOUR MODELLED OBC SHOT GATHER HERE!\n",
    "clip_seis = 1e-7\n",
    "extent_seis = [0.0,xmax/1000,tmax,0.0]\n",
    "\n",
    "plt.subplot(121)\n",
    "plt.imshow(seis_marm_3pt.T, cmap=plt.cm.gray, aspect=2, vmin=-clip_seis, \n",
    "                   vmax=clip_seis, extent=extent_seis)\n",
    "\n",
    "plt.title('3-point operator')\n",
    "plt.xlabel('x [km]')\n",
    "plt.ylabel('t [s]')\n",
    "\n",
    "\n",
    "ax = plt.subplot(122)\n",
    "plt.imshow(seis_marm_5pt.T, cmap=plt.cm.gray, aspect=2, vmin=-clip_seis, \n",
    "                   vmax=clip_seis, extent=extent_seis)\n",
    "ax.set_yticks([]) \n",
    "plt.title('5-point operator')\n",
    "plt.xlabel('x [km]')\n",
    "#plt.ylabel('t [s]')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What we learned:\n",
    "\n",
    "- How to model wave propgation in the complex Marmousi-2 model. Now you can model everything ..."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
